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ABSTRACT 

We present a simple model for the solar differential rotation and meridional circulation based on 
a mean field parameterization of the Reynolds stresses that drive the differential rotation. We in- 
clude the subadiabatic part of the tachocline and show that this, in conjunction with turbulent heat 
conductivity within the convection zone and overshoot region, provides the key physics to break the 
Taylor-Proudman constraint, which dictates differential rotation with contour lines parallel to the axis 
of rotation in case of an isentropic stratification. We show that differential rotation with contour lines 
inclined by 10° — 30° with respect to the axis of rotation is a robust result of the model, which does 
not depend on the details of the Reynolds stress and the assumed viscosity, as long as the Reynolds 
stress transports angular momentum toward the equator. The meridional flow is more sensitive with 
respect to the details of the assumed Reynolds stress, but a flow cell, equatorward at the base of the 
convection zone and poleward in the upper half of the convection zone, is the preferred flow pattern. 
Subject headings: Sun: interior — rotation — Sun: helioseismology 
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1. INTRODUCTION 

Helioseismology has revealed detailed information 
about the internal rotation of the sun. Although the 
radiative interior of the Sun shows a roughly uniform ro- 
tation, the convection zone shows a differential rotation 
with a pole-equator difference of about 30% of the core 
rotation rate. The transition between the two regions is 
confined in a narrow shear layer, the so called tachocline, 
centered at about Q.l Rc^ and havin g a thickness of about 
0.04 i?Q l|Charbonneau et al.iri999|) . Within the convec- 
tion zone the contours of constant angular velocity show 
in mid latitudes an inclination of about 2 5° wit h respect 
to the axis of rotation (Schou et al. 1998i.l2002() . 

There is also robust observational evidence for a 
meridional flow, which is poleward near the surface 
and has a flow velocity around 20 ms^^. This flow 
was first observed through movem ents of acti ve re- 
gions and magnet ic filaments (|Labontc fc Howarc]lll982t 
iTonka et al.lll982l) and later con firmed th rough helioseis- 
mology (see, e.g..'Bra uii fc Faniri998i .Haber et alJl2002l: 
|Zhao & Kosovichcv 200^^ Whereas most of the helio- 
seismic studies focus only on the flow field very close to 
the su rface (using local helioseismology), iBraun fc Fan! 
()1998|) used global helioseismology to extend their anal- 
ysis further down and found a poleward meridional flow 
in the entire upper half of the convection zone. Even 
though the return flow has not been detected yet, it is 
self-evident that it is located in the lower half of the con- 
vection zone because of mass conservation. Since the 
density is significantly larger at the base of the convec- 
tion zone, the fiow amplitude can be very small there, on 
the order of a few ms~^. ^^^^^^^_^^^^^^^ 

Fo llowing initial work bv | Glatzmaier fc Gilmanl (|1982f ) 
and iGilman fc Millen l)1986|) . differential rotation has 
been addressed over the past few decades by full 
spherical shell simulations of compressible convection. 
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Even though recent re sults by iMiesch et al.l l)2000() and 
iBrun fc Toomrd l)2002|) show improvement, the simula- 
tions still have problems in reproducing the radial gra- 
dient of the differential within the convection zone. Al- 
though the models predict the magnitude of the differ- 
ential rotation (latitudinal variation) correctly, the con- 
tour lines of constant angular velocity are still very close 
to the Taylor-Proudman state with isolines parallel to 
the axis of rotation. The three-dimensional simulations 
also have problems in generating the observed meridional 
flow, which shows in the upper half of the convection zone 
a poleward flow of an amplitude of about 20ms~^. In 
contrast to observations, the three-dimensional simula- 
tions typically show several smaller flow cells of opposite 
sign and large temporal variability. 

A completely different approach is based on ax- 
isymmetric mean field models that parametrize the 
turbulent angular m o mentu m transport (A-effect; 
[Kitchatinov fc Riid iged I1993J) . As discussed by 

iRudigcr et al. (199^, these models typically show solu- 
tions close to the Taylor-Proudman state unless a very 
large value of the turbulent viscosity is used. 

[Kitchatinov fc Riidigcr ( 1995) showed that an 
anisotropic convective energy transport can produce a 
latitudinal variation in the temperature that is large 
enough to overcome the Taylor-Proudman constraint 
and allow for solar-like differential rotation while assum- 
ing reasonable values for the turbulent viscosity. The 
anisotropy of the energy flux results from the rotational 
influence on the turbulence mo difying the turbulent 
heat diffusivities IjKitchatinov et al. 19 9#. Recen tly 
this model has been applied bv lKiiker fc Stixl l)200H) to 
study the evolution of the solar differential rotation of 
the Sun from the pre-main-sequence Sun to the present 
Sun by solving the momentum equations together with 
an mixing-length approach for the convective energy 
flux. The role of latitudinal variations of the entropy 
for the diffe rential rotation was also studied bv iD^irnevI 
(119991 120031) within the framework of mean field models. 
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In this paper we present a model that is along the lines 
of the axisynimetric mean field n iodels mentioned above. 
Whereas the investigations of iKitchatinov fc Riidieeil 
Ijl995(l and Kiikcr & Stix (2001) focused on the role of 
an anisotropic convective energy flux, we focus our in- 
vestigation on the role of a subadiabatic tachocline for 
the Taylor-Proudman balance of the differential rotation. 
The aim of this paper is not to present a complete model 
for the solar differential rotation, but rather a simple 
approach to investigate the specific question of how a 
subadiabatic tachocline, in conjunction with turbulent 
heat conductivity within the convection zone and over- 
shoot region, can break the Taylor-Proudman constraint 
which requires a differential rotation constant on cylin- 
ders in case of an isentropic stratification. 

2. MODEL 

In this investigation we use a simplified model for 
studying differential rotation and meridional circulation 
in the solar convection zone. The basic assumptions un- 
derlying this approach are as follows. 

1. Axisymmetry and a spherically symmetric refer- 
ence state. 

2. All processes on the convective scale are param- 
eterized, leading to turbulent viscosity, turbulent 
heat conductivity, and turbulent angular momen- 
tum transport. 

3. The equations can be linearized assuming gi <^ go 
and pi <^ po, and the reference state is assumed 
to be spherically symmetric. Here go and po de- 
note the reference state values, whereas gi and pi 
are the perturbations caused by the presence of dif- 
ferential rotation. The equations we solve are the 
fully compressible, linearized, axisymmetric hydro- 
dynamic equations. 

4. The entropy equation includes only perturbations 
associated with differential rotation. Therefore, the 
reference state is assumed to be in an energy flux 
balance. It is further assumed that the effect of 
convection (in the convection zone and the over- 
shoot region) on entropy perturbations associated 
with the differential rotation is purely diffusive. 

5. The tachocline at the base of the solar convection 
zone is forced by a uniform rotation boundary con- 
dition at r = 0.65 i?0. 

We emphasize that this model is not intended to be a 
complete solar convection zone model, since fundamental 
processes required for differential rotation such as turbu- 
lent angular momentum transport are parameterized. 

This model is also intended to be a basis for a 'dy- 
namic' flux-transport dynamo including the j X B force 
feedback on meridional flow and differential rotation (fu- 
ture work). 

2.1. Basic equations 

For this model we use the axisymmetric, fully com- 
pressible hydrodynamic equations. Since the perturba- 
tion of pressure and density caused by the differential 
rotation are small compared with the reference state val- 
ues [gi/go ~Pi/po ~ (Ar2i?0/cs)^ ^ 10"^], we linearize 



the equations, assuming gi ^ go and pi ^po- Since we 
do not use the anelastic approximation here (see section 
12.61 for more details), we keep the time derivative in the 
continuity equation: 
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We use here the dimensionless entropy s = hi(ji/g'^), 
meaning that the entropy equation Eq. |SJ| was made 
dimensionless by division through c„ = (7 — l)~^i?//i. 

The third term on the right-hand side of Eq. Q 
describes the effects of a non-adiabatic reference state, 
where (5 = V — Vad is related to the gradient of the ref- 
erence state entropy through: 

dso jS 
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The fourth term in Eq. Q considers the energy trans- 
fer through Reynolds stress that will be defined in de- 
tail in the following paragraph. The last term describes 
turbulent diffusion of entropy perturbation within the 
convection zone, where nt denotes the turbulent thermal 
diffusivity. We neglect in Eq. Q the contribution of the 
radiative energy ffux, since even for overshoot values the 
turbulent heat conductivity exceeds the radiative one by 
several orders of magnitude. 
The viscous force F follows from 
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with the Reynolds stress tensor 



Rik = -go < v^v^. >= vtQo ( Eik - -5ik divt) + Ajfc 

(12) 
Here Eik ~ Vi-k + Vk-i denotes the deformation tensor, 
which is given in spherical coordinates by 
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whereas Aik denotes the non diffusive Reynolds stresses, 
which are responsible for driving differential rotation. 
We will discuss these terms later. 

The amount of energy that is converted by the 
Reynolds stress is given by: 



Q — 2_^ 1^ Eik Rik ■ 



(19) 
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Q contains a heating term resulting from the dissipation 
of kinetic energy through the dissipative contribution to 
the Reynolds stress (terms proportional to Eik) and a 
cooling term resulting from the energy transfer intro- 
duced by the nondiffusive transport term proportional 
to Aik responsible for maintaining the differential rota- 
tion. The latter is in general the dominant term. The im- 
portance of the term Q for a stationary solution becomes 
more apparent if we transform the entropy equation to an 
equation of the quantity qoTqSi, which better represents 
the energy perturbation associated with the entropy per- 
turbation. In the case of a stationary solution, we have 
div(go'u) — and we can rewrite the entropy equation 
Eq. (O in the form (assuming \S\ = |V — Vad| ^ 1) 



div {v goToSi - KtQoTo gradsi) 
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The terms that appear in this equation are the diver- 
gence of the energy flux (left-hand side), a source term 
that considers the viscous heating and the buoyancy work 
(first term, right-hand side), and the source term that 
arises from the nonadiabatic stratification (second term, 
right-hand side). The first term redistributes goToSi but 
does not provide a net source, since the flux across the 
boundaries vanishes with the boundary conditions we 
use. The same applies to the last term in the case of 
a stationary solution, since the horizontal mean of the 
radial mass flux Vrgo has to vanish and S and g show no 
latitudinal dependence. The only net source is the second 
term, which contains the Reynolds stress and buoyancy 
work. In the case of a stationary solution the volume 
integral of this term has to vanish, which means that the 
energy extracted through the A-effect from the reservoir 



of internal energy (through feedback on convective mo- 
tions) returns through viscous heating and work of the 
meridional flow against the buoyant force. 

2.2. Background stratification 
For the background go, Pa, and Tq we use an adiabatic 



hydrostatic stratification assuming an ~ r 
of the gravitational acceleration given by 
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where Tbc, Pbc, and ghc denote the values of temperature, 
pressure and density at the base of the convection zone 
r = Tbc. Here iJbc = Phc/{ghcghc) is the pressure scale 
height and g^c the value of the gravity at rbc- In the 
following we use rbc = 0.71 Rq pbc = 6 x 10^^ Pa, gbc — 
200kgm-3, .gbc = 520 ms'^, and Rq = 7x lO^m, which 
results in H^c = 0.0825 i?0. 

2.3. Superadiabaticity profile 

For the superadiabaticity S we assume the following 
profile: 

5 = Sconv + 7r(<5os - ^conv) I 1 " tanh ( 

^ V V "trail 



where 



(^coiiv = (5topexp(- 
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Here Stop, (^cz, and Sos denote the values of supera- 
diabaticity at the top of the domain r = rmax, in the 
bulk of the convection zone, and in the overshoot re- 
gion, respectively. In addition rguh denotes the radius at 
which the stratification within the convection zone turns 
weakly subadiabatic (because of nonlocal effects), and 
''tran dcuotcs the radius of transition towards stronger 
subadiabatic stratification (the overshoot region). The 
parameters dtop and dtran determine the steepness of the 
transition toward large superadiabaticities at the top of 
the domain and towards the overshoot region, respec- 
tively. Fig. n shows the profiles of the superadiabaticity 
we use later in our models. 

Nonlocal mixing-length models show a transition from 
subadiabatic to superadiabatic values typically between 
r = 0.75 and 0.8 Rq depending on the assumed mixing- 
length parameter (PidatcUa & Stix 1986; Skalcv & Sti3 
Il991|) . For larger degrees of non- locality the subadia- 
batic fraction of the convect ion zone could be even larger 
l|SDruitlll997t iRempell 12001 . Mixing-length models pre- 
dict a variation of S within the convection zone by several 
orders of magnitude; however, most of this variation oc- 
curs very close to the surface layers, which are difficult 
to resolve in a mean field model. Since our model cap- 
tures only the large-scale fiows, we also have to make 
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Fig. 1. — Profile of superadiabaticity used in tlie models with 
adiabatic convection zone (solid line), model 8 (dotted line), model 
9 (dashed line), and model 10 (dashed-dotted line). In all cases we 
have (5os = -1.5 X IQ-^, rtran = 0-725 Rq, and dtran = 0.0125 Rq. 
Common parameters of the models 8 to 10 are dtop = 0.0125 i?©, 



5;,^ = 3 X 10" 

'sub 



and <5to 



: 3 X 10" 



Case 8 shows a profile with 
= 0.8 ij0, case 9 with r■g^^, = 0.825 Aq, and case 10 with 
^sub = 0.85 iJ©. Shown on the vertical axis is |5| on a logarithmic 
scale. The singularities at r = 0.8, 0.825, and 0.85 -Rq indicate 
where 5 changes sign in the models with a nonadiabatic convection 



sure that the Rayleigh number stays sub critical for con- 
vection within the domain for the thermal conductivity 
and viscosity we use. 

In this investigation we keep the following parameters 
fixed: dtop = ^tran = 0.0125 i?0 and rtran = 0.725 i?©. 
Choosing a value of (5top = 3 x 10^^ we have to use a 
thermal diffusivity of around 5 x 10^ m^ s~^. Much larger 
values of b would require unreasonably large values of 



the turbulent heat conductivity. Below r = rtran we as- 
sume overshoot type values of the superadiabaticity, in 
the range — 10~^ to —10""'. We extend our domain down 
to r = 0.65 i?o, which also includes the radiative interior 
with values of 5 ~ —0.1. Since our diffusivity profile 
(see next subsection) drops significantly below rbc, the 
entropy perturbation generated at the lower boundary 
does not influence the values within the convection zone 
and has therefore no influence on the differential rota- 
tion profile. More important for this model is the overlap 
between the subadiabatic region and the thermal diffu- 
sivity profile, since this determines to what extent an 
entropy perturbation can spread from the subadiabatic 
region into the convection zone. Since the representation 
of a strongly subadiabatic layer would cause numerical 
difficulties (see section IT^ for further details) we decided 
to exclude this layer and use overshoot values for b down 
to 0.65 Rq. The effective thickness of the overshoot is de- 
termined by the overlap with the thermal conductivity 
profile, which we define in the next subsection. Since the 
heat conductivity is basically zero below rbc = 0.71 i?0, 
the effective thickness of the overshoot region is about 
2 (itran ~ 20 Mm, which is in the range of predictions 
of nonlocal mixing- l ength models (PidatcUa & Stix 1986; 
|Skalev& Stix 1991). 

Since all these values are small compared with Vad, we 
still can use the adiabatic reference state for pressure, 



density, and temperature as given above in Eqs. H21|l to 
(ESI)- 

2.4. Diffusivity profiles 

For the turbulent viscosity and thermal conductivity 
we assume constant values within the convection zone 
and the radiative zone with a transition smoothed by a 
hyperbolic tangent function. We assume that the diffu- 
sivities only depend on the radial coordinate: 
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where vq and kq denote the values of the turbulent dif- 
fusivities within the convection zone and a^i/ specifies 
the values of the turbulent diffusivities at r = rtran 
(ofKi/i'o, aKi^«^o)i where the transition to the overshoot 
region takes place in our model. The profile function 
/c ensures that the diffusivities drop significantly toward 
the radiative interior at rbe- For the width of this transi- 
tion we use dbe = 0.0125 i?Q. Since both diffusivities are 
of turbulent origin, we use the same radial profile. 

In this model we do not include a sophisticated the- 
ory for the tachocline but force a tachocline through a 
uniform rotation boundary at r = 0.65 i?0. As a con- 
sequence the tachocline is a viscous shear layer between 
the convection zone and the lower boundary condition. 
We therefore have to maintain a sufficient amount of vis- 
cosity in the radiative interior to allow for the formation 
of this shear layer in a reasonable amount of time. In the 
following we use the profile defined by Eq. (|27|l for the 
computation of the turbulent angular momentum trans- 
port (last term in Eq. ^21) but set the viscosity used 
for the diffusive terms of the Reynolds stress (first two 
terms in Eq. J21) to 2% of the convection zone values. 
For the heat conductivity we use in the radiative interior 
0.2% of the convection zone values. 

2.5. Parameterization of turbulent angular momentum 

transport 

The terms relevant for the differential rotation are: 






-L{r,0) cos(6i 
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(31) 
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where L(r, 9) denotes the amplitude of the angular mo- 
mentum flux, whereas A(r, 9) describes the inclination of 
the flux vector with respect to the axis of rotation. 

Symmetry considerations require that A^^ is symmet- 
ric and Ag^ is antisymmetric across the equator. To sat- 
isfy this constraint, A(r, 9) and L(r, 9) need to be anti- 
symmetric functions with respect to the equator. Since 
we solve our model only in the northern hemisphere, we 
specify in the following discussion always values for A 
and L in the northern hemisphere. To maintain the 
proper symmetry for a full-sphere simulation, values for 
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the southern hemisphere would have to be chosen by re- 
flection across the equator. 

The setting L > 0, A = corresponds to a flux directed 
downward and parallel to the axis of rotation; L > 0, 
A — ~9 to a radially inward flux, and L > 0, A = 7r/2 — 
to an equatorward flux. The setting L r^ sin 6 cos 9 and 
A = rec overs the limit of fast rotation found for the 
A-effect bv lKitchatinov fc Riidigal 1)19951) . 

In the following discussion we will for the amplitude of 
the angular momentum flux the expression 



/(r,6l) = (sin6')"cos6'tanh 
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The full radial dependence of the angular momentum 
flux is obtained by multiplication by vtQo- The angular 
momentum flux drops below rtran because of the drop in 
turbulent viscosity. In most of the following discussion 
we require a vanishing angular momentum flux at the 
top boundary as expressed in Eq. (|33|) . where we use a 
value oi d — 0.025 Rq for the transition layer at the top. 
The exponent n determines the latitude at which the 
flux peaks. For A > 0, the value of n needs to be larger 
than 2 to ensure the regularity of the divergence of the 
Reynolds stress close to the pole. For the direction of 
the flux determined by A we discuss two distinct cases: 
A = 15° (transport nearly parallel to axis of rotation) 
and A = 90° — 9 (transport in latitude only) in order to 
evaluate the sensitivity of the model with respect to this 
parameterization. 

2.6. Numerical procedure 

We are interested here in the stationary solution of 
Eqs. Q to |SJ for a given parameterization of the tur- 
bulent angular momentum transport Eq. 1)31)) . (|32() . A 
very natural way to relax the system is to use the tempo- 
ral evolution; however, because of the low Mach number 
of the expected flows, a direct compressible simulation is 
problematic. For an expected meridional flow velocity of 
a few ms~^ the Mach number is around 10~^ (because 
of axisymmetry, the much faster differential rotation flow 
does not enter the time step limit). Without leaving 
the regime of highly subsonic flows and therefore with- 
out changing the physical properties of the solution, it is 
possible to speed up the relaxation process signiflcantly 
by increasing the base rotation rate f2o- 

Using the following transformation for the independent 
parameters of the equations. 
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and the following transformation for variables. 
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Eqs. © " © remain unchanged, meaning that if 
{gi, Vr, Vg, rii, si} is a solution for the parameters 



{ilo, vt, Kt, 6} then {^gi, Ci'r, C^e, C^i, C^si} is a 
solution for the parameters {C^o, C^t) C'*t7 C^^}- How- 
ever this transformation changes the equation of conti- 
nuity to: 

dgi 1 



— -f — div(go -y) - , 



(35) 



meaning that the time evolution is changed, but the sta- 
tionary solution remains unchanged. The pre factor C~^ 
in the equation of continuity corresponds to the increase 
of the Mach number of the flow by a factor of C- In the 
following we use a value of C = 100, which corresponds to 
an increase of the Mach number of the meridional flow 
in the bulk of the convection zone from 10~^ to 10"'^. 
Therefore even the time evolution is only marginally af- 
fected, since the solution stays in the regime of highly 
subsonic flows. We computed solutions with different 
values of ^ in order to quantify the influence of this trans- 
formation. We found differences between a solution com- 
puted with C = 10 and 100 on the order of a few percent 
(mainly in the magnitude of the differential rotation close 
to the pole). 

This approach is very similar to simulations of rising 
magnetic flux tubes in the solar convection zone, which 
were made fully compressible by assuming a value oi (3 = 
Pgas/Pmag On the Order of 100 instead of 10^ to overcome 
a severe time step constraint. As long as the relevant flow 
velocities remain sufficiently subsonic, the results are not 
affected significantly. 

The only drawback of this approach is that it is im- 
possible to represent large values of the subadiabaticity 
as found in the radiative interior of the Sun. Because 
of the scaling of S with C^, a value oi S = —0.1 would 
correspond to S — —1000 in a solution with a base ro- 
tation increased by a factor of 100, which is physically 
impossible. However, overshoot-like values of 6 ^ — 10~^ 
can be treated without any problem. As we explained 
earlier in subsection l2.3l and 12 . 41 the radiative interior is 
not of great importance, since entropy perturbations cre- 
ated there cannot influence the convection zone we are 
primarily interested in. 

Since our model also includes a signiflcant time step 
constraint because of the large turbulent diffusivities 
in the convection zone, an anelastic approach would 
not provide much advantage unless all diffusivities were 
treated implicitly. We therefore decided to solve the 
equations with a faster explicit scheme using the pro- 
cedure outlined above. We tested an anelastic version 
of the code and found convergence problems of the pres- 
sure solver related to the uniform rotation boundary that 
we impose at r = 0.65 i?©. For a different choice of 
boundary conditions we found very good agreement be- 
tween the anelastic version and a solution computed with 
C = 100 as described above. 

We solve Eqs. O - © with a MacCormack scheme us- 
ing alternating upwind and downwind differencing, which 
is second order in space and time. The computational 
domain extends in latitude from equator to pole and in 
radius from r = 0.65 to 0.985 i?©. We use the appropri- 
ate symmetry boundary conditions at equator and pole 
and closed boundaries in radius. At the bottom bound- 
ary (r = 0.65 Rq) we enforce a uniform rotation; the top 
boundary is stress-free for the angular velocity (i?^^ — 0). 
At both radial boundaries we use stress-free boundary 
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TABLE 1 

Significant parameters of simplified model 



(f) component of the vorticity under the assumption of 
small deviations from adiabaticity (|V — Vad| ^ 1) 



case 


A 


n 


dKU 


CIku 


l5os 


^sub 


1 


15° 


2 


0.025 


0.1 


-1.5 X 10-5 




2 


90° -e 


2 


0.025 


0.1 


-1.5 X 10-5 




3 


15° 


2 


0.025 


0.1 







4 


15° 


4 


0.025 


0.1 


-1.5 X 10-5 




5 


15° 


2 


0.025 


0.1 


-3 X 10-5 




6 


15° 


2 


0.05 


0.1 


-1.5 X 10-5 




7 


15° 


2 


0.025 


0.025 


-1.5 X 10-5 




8 


15° 


2 


0.025 


0.1 


-1.5 X 10-5 


0.8 


9 


15° 


2 


0.025 


0.1 


-1.5 X 10-5 


0.825 


10 


15° 


2 


0.025 


0.1 


-1.5 X 10-5 


0.85 



Note. — Cases 1-7 have adiabatic convection zones and 
cases 8-10 have nonadiabatic convection zones with, the 
common parameters (5cz = 3 X 10"'' and (Stop = 3 
but different values of r^^\,, as shown in the table. 



X 10" 



conditions for the velocity and set the derivative of si to 
zero. Since we focus in this study on the large-scale flow 
fields, a moderate resolution of around 108 grid points 
in radius and 72 grid points in latitude is sufficient. We 
tested our cod e by reproducing the result presented by 
iRtidJger et all 1)19981 their Fig. 1 ). 

3. RESULTS 

Since our simplified model contains parameterizations 
of crucial processes, we have to evaluate carefully the 
dependence on particular choices of these parameters. 
In Tabled we have summarized the parameters that we 
discuss in this section. There are additional model pa- 
rameters, which do not have a significant influence on 
the solution. These are rtram Aq, Vq, and kq- Here rtran 
specifies where the superadiabaticity turns from convec- 
tion zone to overshoot values. In the following we use 
'"tran — 0.725 i?©, which is a reasonable choice for the 
Sun. As long as rtran > ^bc + dhc the influence on the 
solution is small, since it is more the profile of vt and Kt 
in relation to rtian that matters. We have therefore in- 
troduced in Eqs. H27|l and (|28f) the parameter a^iy, which 
specifies the profile relative to rtran- Similarly Aq spec- 
ifies the magnitude of the nondiffusive Reynolds-stress 
(A-effect). Typical values for Aq are on the order of unity. 
Except for case 2, in which we use Aq = 0.4, in the fol- 
lowing we use a value of Aq = 0.8. We have chosen Aq 
such that the magnitude of the differential rotation is 
close to solar-like. 

For the diffusivities vq and kq we assume in the follow- 
ing discussion i^g = ^o = 5 x lO^m^s"^. For cases with 
a superadiabatic convection zone we have to increase the 
value of Kg to 5 X 10^ m^ s"^ above rgub in order to avoid 
convective instability. 

3.1. General solution properties: differential rotation 

The key ingredient in this model is the inclusion of 
a subadiabatic tachocline beneath the convection zone, 
which is enforced in this model through the uniform ro- 
tation lower boundary condition. Within the subadia- 
batic region (r < 0.725 Rq) the differential rotation is 
balanced by a latitudinal entropy gradient. Taking the 
curl of the meridional momentum equation yields for the 






r sm f 



9 dsi 

jr do 



'^--1^ (36) 



with Q = Qq + fli] the bracket denotes viscous terms 
and vorticity transport terms, which are not important 
for the following discussion. 

Starting initially with si =^ 0, the turbulent angu- 
lar momentum transport leads to a negative values of 
dVP' /dz in high latitudes, which enforces a negative value 
of (jj^. A negative value of w^ corresponds to a counter- 
clockwise meridional flow in the tachocline, which shows 
a negative radial velocity at high latitudes and a positive 
velocity at low latitudes. Because of the subadiabatic 
stratification, this results in a positive entropy pertur- 
bation in high latitudes and a negative entropy pertur- 
bation in low latitudes, as shown in Fig. Ob). Since 
the resulting negative value of dsi/dO can compensate 
for the also negative value of dVi} jdz, an equilibrium is 
reached finally. An additional source for the entropy per- 
turbations comes from the meridional flow driven in the 
convection zone and penetrating to some extent into the 
subadiabatic overshoot region. For the parameters used 
in this model both effects are of roughly the same order 
of magnitude. In our model most of the tachocline shear 
is located below the base of the convection zone, whereas 
helioseismic inversions find more o verlap between the 
tachocline and the convection zone IjCharbonneau et alJ 
I1999|) . Our model has therefore most probably the ten- 
dency to underestimate the entropy perturbation in the 
overshoot region caused by the value of d^ jdz. 

Because of the turbulent thermal heat conductivity 
this entropy perturbation can spread into the convec- 
tion zone and therefore also balance there a differen- 
tial rotation that deviates from the Taylor-Proudman 
state with J7-contours parallel to the axis of rotation. 
We want to emphasize that the total entropy sq + si in 
the overshoot region is still smaller than in the convec- 
tion zone. The physical reason for this spread is that 
the convection tries to maintain the same radial entropy 
gradient at all latitudes (if we do not consider possible 
rotational anisotropy). Since the overshoot region pro- 
vides the entropy boundary condition for the convection 
zone, a latitudinal variation of entropy in the overshoot 
region is tr ansported by convection into the convection 
zone. iStixl l)1981|) computed response functions for the 
temperature, velocity, and flux perturbations within in 
the framework of the mixing-length approach and found 
that the screening effect of temperature perturbations 
is very weak, meaning that temperature (entropy) per- 
turbations at the base of the convection zone should be 
transmitted through the entire convection zone. 

The magnitude of the entropy perturbation in the con- 
vection zone depends on the overlap between the thermal 
conductivity profile and the subadiabaticity profile. For 
most models we use a parameter of a^iy = 0.1, which 
means that the thermal diffusivity drops to 10% of its 
convection zone values at r — rtran = 0.725 i?©, but 
smaller values (ckki/ — 0.025) also work if the value of 
Sos is slightly increased (see cases 5 and 7). The ther- 
mal heat diffusivity in a nonlocal mixing-length model, 
^ VconvHp, would lead to larger values within the over- 
shoot region because of the large value of the pressure 
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Fig. 2. — Contours of (a) differential rotation and (b) entropy perturbation for case 1. Solid lines indicate positive values. The entropy 
perturbation that originates in the subadiabatic tachocline and spreads because of thermal conductivity into the convection zone prevents 
the Taylor-Proudman state (contours parallel to axis of rotation) for the differential rotation from developing, (c) Contours of differential 
rotation for case 3. This case is similar to case 1, except that si = 0. As a consequence the contour lines of constant fl are aligned with 
the axis of rotation. 



scale height and still significant velocities ~ 10 ms"^ in 
the overshoot region. 

The temperature perturbation associated with the en- 
tropy perturbation shows a magnitude of about 5K 
throughout the convection zone for most of our mod- 
els. Whereas the entropy perturbation drops monoton- 
ically from pole to equator, the temperature perturba- 
tion close to the surface reaches a minimum around mid 
latitudes followed by a slight increase toward the equa- 
tor. At the base of the convection zone the temperature 
also shows a monotonic decrease from the pole toward 
the equator. The different behavior of entropy and tem- 
perature is due to the pressure perturbation within the 
convection zone, which also contains an adiabatic con- 
tribution (first term on right-ha nd side of Eg. E|). A 
similar pattern was also found bv lBrun fc Toomrd (|20fl2fl 
in three-dimensional simulations; however, the physical 
reason might be different. 

In Fig. |2] we show the contours of fl and the re- 
lated entropy perturbation si for case 1. Since the en- 
tropy perturbation is concentrated in higher latitudes 
(where dil^/dz also peaks), the deviations from the 
Taylor-Proudman state are largest in high latitudes. 
Whereas the differential rotation shows close to the pole 
J7-contours perpendicular to the axis of rotation, the Q- 
contours are more aligned with the axis of rotation close 
to the equator. 

In Fig. El c) we show for reference purposes a solu- 
tion with si = but, apart from that, exactly the same 
parameters as case 1. Without the effect of the subadi- 
abatic tachocline the solution shows differential rotation 
with cylindrical Q contours, even though we still impose 
the uniform rotation boundary condition at r = 0.65 Rq. 

3.2. General solution properties: meridional flow 

Since the assumed turbulent angular momentum trans- 
port has in general the tendency to drive a differential 
rotation that differs from the profile that would balance 
the right-hand side of Eq. 1)36(1 , the differential rotation 
shows a small perturbation around this state. As a con- 
sequence, the Coriolis force related to this perturbation 
cannot be precisely balanced by a combination of pres- 



sure gradient and buoyancy and therefore drives a merid- 
ional flow (examples of this flow can be seen in Fig. Q . 
This flow grows until the additional angular momentum 
transport sufficiently limits the perturbation of f2 so that 
the remaining weak unbalanced Coriolis force can be bal- 
anced by viscous stress associated with the meridional 
flow. For the parameterization of the Reynolds stress 
used in our model we typically get a counterclockwise 
flow cell (equatorward at base and poleward at surface) . 
This is strongly related to the radial component of the 
turbulent angular momentum flux, which is assumed to 
be directed inward in most of our models. In the absence 
of a meridional flow, it would increase the rotation rate 
at the base and decrease the rotation rate at the surface, 
which leads to counterclockwise Coriolis forces. 

The meridional flow typically closes above r = 0.71 Rq 
for two reasons. (1) The subadiabatic stratiflcation sup- 
presses radial motions. (2) The Reynolds stress (~ vt) 
that indirectly drives a meridional flow through a change 
of il drops significantly below r = 0.71 i?©. The merid- 
ional fiow will never vanish completely in the subadi- 
abatic region, since some radial motion is required to 
maintain the entropy perturbation against diffusive de- 
cay. However, the flow velocities are very small com- 
pared with the values within the convection zone. For 
our choice of a value of^^— 10~^, the value of ve is on 
the order of a few cms~^. For values of ^ ~ — 10~^ as 
found in the radiative interior, these values would drop 
by at least 3 orders of magnitude. The limits for the 
penetration of the meridional flow below the base of the 
convection zone ar e even more stringent in this model 
than the analysis of lGilman fc MieschI l|2004|) showed for 
overshoot-like values of 6. The main reason for this is 
the inclusion of the feedback of the meridional flow on 
the differential rotation through the transport of angu- 
lar momentum. Any signiflcant equatorward flow in a 
region with low turbulent diffusivity would cause a ret- 
rograde zonal ffow that would oppose the flow through 
the Coriolis force. Beside the effect of subadiabaticity 
to suppress radial motions, the angular momentum con- 
servation in a low-diffusivity region yields an additional 
effect suppressing latitudinal motions. Only if there is a 
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Fig. 3. — Influence of different parameterizations of the turbulent angular momentum transport on differential rotation, (a - c) Contours 
of the differential rotation, where solid lines indicate regions rotating faster than the core, (d - f) Profile of the differential rotation as a 
function of radius for the latitudes 0°, 15°, 30°, 45°, 60°, and 90°. Differential rotation is shown for (a, d) case 1 (A = 15°), (b, e) case 2 
(A = 90° — d), and (c, f) case 4 (A = 15° and n = 4). The profile of the differential rotation is rather insensitive to changes in the Reynolds 
stress; however, the amplitude changes, especially in high latitudes. 



process that disturbs an equilibrium solution of Eq. H36|l 
can there be a significant meridional flow below the base 
of the convection zone. 

3.3. Dependence on parameterization of angular 
momentum transport 

Figs. O and 21 show the sensitivity of the solution with 
respect to different parameterizations of the turbulent 
angular momentum transport. In cases 1 and 2 the am- 
plitude of the angular momentum flux has the same pro- 
file in radius and latitude; however, the direction of the 
flow is changed. In case 1 the angular momentum is 
transported almost parallel to the axis of rotation with 
a 15° inclination angle; in case 2 the angular momen- 
tum flux has only an equatorward latitudinal component 
(A = 90° —9). In both cases the amplitude of the an- 
gular momentum flux, Aq, was adjusted such that the 
equatorial rotation rate is the same (Aq = 0.8 in case 1 
and Aq = 0.4 in case 2). The variation in the required 
amplitude follows from the fact that an angular momen- 
tum transport perpendicular to the axis of rotation is 
the most efficient way to speed up the rotation at the 
equator, whereas a transport parallel to the axis of ro- 
tation would have no effect at all. A comparison of Fig. 
13 panel a) and b), shows that the profile of the differ- 



ential rotation is not very sensitive to the change in the 
direction of the angular momentum. 

Fig. O panels c) and f) show case 4, which is simi- 
lar to case 1, but uses a turbulent angular momentum 
transport confined closer to the equator (n = 4 instead 
of n = 2). Whereas the profile of the differential rota- 
tion is only marginally changed, the magnitude of the 
differential rotation in high latitudes is reduced. 

Unlike fi, the meridional flow is much more sensitive 
to the details of the angular momentum transport. Fig. 
01 panels a) to c), show the streamlines of the meridional 
flow, where solid lines indicate counterclockwise flows 
(poleward at the surface and equatorward at the base 
of the convection zone). Fig. 01 panels d) to f), show the 
latitudinal component of the meridional flow velocity at 
45° latitude. All cases show a dominant counterclockwise 
cell; however, each case shows different flow amplitudes. 
The cases 1 and 3 with a significant angular momentum 
flux along the axis of rotation show flow velocities by 
a factor of around 3 larger, which is in part caused by 
the larger value of Aq required in these cases to obtain 
the same equatorial rotation rate. Case 2 with angular 
momentum flux in latitude only also shows a weak re- 
verse circulation cell above 60° latitude. The second cell 
is driven by buoyancy resulting from the higher value of 
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Fig. 4. — Influence of different parameterizations of the turbulent angular momentum transport on meridional flow, (a - c) Stream 
function, where solid lines indicate counterclockwise flows, (d - f) Radial profile of Vg at 45° latitude. Meridional fiow is shown for (a, 
d) case 1, (b, e) case 2, and (c, f) case 4. Unlike the differential rotation, the meridional flow is very sensitive to the direction of the 
angular momentum transport (compare cases 1 and 2). In case 4 the confinement of the Reynolds stress to lower latitudes also confines 
the meridional fiow to lower latitudes. 



the entropy close to the pole. This second cell is not 
visible in case 1, since the radial component of the an- 
gular momentum flux has a strong tendency to drive a 
counterclockwise meridional flow. The radially inward 
flux of angular momentum at high latitudes (because of 
the assumed 15° inclination with respect to the axis of 
rotation) increases the value of dQ'^/dz in high latitudes, 
which leads to a counterclockwise meridional flow accord- 
ing to Eq. H3t)|l . There is also a very weak indication of 
this second cell in case 4, in which the larger value of n 
confines the meridional flow to lower latitudes. 

The main reason for the different sensitivities of 
meridional flow and differential rotation with respect to 
changes in the Reynolds stress follows from the fact that 
the differential rotation is mainly determined by the bal- 
ance expressed in Eq. H36|l . Since the entropy pertur- 
bation is not directly affected by a changing Reynolds 
stress, the profile of fl also changes little. The merid- 
ional flow, on the other hand, results from an imbalance 
of Eq. (|36() and is therefore the result of a small difference 
between large forces. As consequence, the sensitivity of 
the meridional flow to changes in the Reynolds stress is 
much larger. 

We focus here on two distinct cases for the value of 
A. We want to mention that any choice for A > 0° and 



< 90° leads to differential rotation close to case 1 (large 
values of A are closer to the Taylor-Proudman state in low 
latitudes); however, the meridional flow changes signifi- 
cantly. Changing the value of A requires an adjustment 
of Aq ~ 1/ sin A in order to keep the magnitude of the dif- 
ferential rotation fixed. As a consequence, models with 
smaller values of A have larger meridional flow velocities 
than models with larger values of A. Using a value of 
A > 45° also leads to a more complicated flow structure, 
which shows a clockwise flow pattern within most of the 
convection zone. Using values of A < 15° leads to so- 
lutions very similar to case 01, but with a significantly 
larger amplitude of the meridional flow velocity. 

3.4. Dependence on subadiabaticity, viscosity, and 
thermal conductivity profile 

Fig. |31 compares models with identical parameteriza- 
tion of the Reynolds stress but different parameters for 
the subadiabaticity and the profiles of viscosity and heat 
conductivity. Panels a) and d) show how different val- 
ues of the subadiabaticity affect differential rotation and 
meridional flow. In both panels the case 1 is shown as a 
reference (dashed line). Case 5 with a value of (5bc twice 
as large shows an increase of the differential rotation by 
about 20%. Whereas the differential rotation remains 
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Fig. 5. — Influence of subadiabaticity, viscosity, and conductivity on solution, (a, d) Case 5 with a subadiabaticity increased by a factor 
of 2, leading to larger equator-pole difference in Q but nearly no change in the meridional flow pattern, (b, e) Case 6 with an increased 
value of d^i/ , leading to significant reduction in the meridional flow speed, (c, f ) Case 7 with a decrease in a^v Decreasing the overlap 
between the diffusivity profile and the subadiabaticity profile is very similar to decreasing the value of <5, except that the meridional return 
flow at the base of the convection zone is located at a different depth. In all panels we have shown case 1 as a dashed line for reference. 



unchanged at the equator, the polar values of O decrease 
significantly. On the other hand, the meridional flow is 
only marginally affected. A larger value of She leads to a 
larger entropy perturbation, which can balance a larger 
differential rotation because of Eq. H36|l without requir- 
ing a change of the meridional flow. A similar result can 
be obtained by lowering kq and keeping She constant. 

In panel b) and e) we show a solution with a factor of 
2 larger width of the transition in the diffusivity profile 
duv (case 6). The most striking change is visible in the 
meridional flow pattern, where the return flow speed at 
the base of the convection zone is reduced by nearly a 
factor of 2. The amplitude of the differential rotation is 
only marginally affected. 

Panel c) and f) shows case 7 with the parameter a^iy 
decreased by a factor of 4 (value of the diffusivities at 
J'tran)- The meridional flow pattern shows a reduction 
of the penetration into the subadiabatic layer because of 
the change in the Reynolds stress that drives this flow. 
The magnitude of the differential rotation at high lati- 
tudes decreases in response to the reduced spread of the 
entropy perturbation into the convection zone. Combin- 
ing a lower value of a^jy with a larger value of the sub- 
adiabaticity in the overshoot region as shown in panel 
(a) and c) would compensate for this effect and provide 



nearly the same solution for the differential rotation (the 
penetration depth of the meridional flow would be still 
different). 

So far we have discussed the influence of different pro- 
files of viscosity and thermal diffusivity, but the mag- 
nitudes set by vq and kq have been left constant. We 
mentioned above that a change in kq has the opposite 
effect as a change in ^bc- Whereas an increase in She 
increases the entropy perturbation and through the bal- 
ance Eq. (|36|l the differential rotation, an increase in 
Kq decreases both entropy perturbation and differential 
rotation. An increase in vq increases the magnitude of 
the meridional flow and, through the term r^ Sv^ in Eq. 
10, also the entropy perturbation, which also results in 
an increase in differential rotation. However, if vq and kq 
are changed together, the latter effect is compensated for 
by the larger thermal diffusivity and the differential ro- 
tation stays the same, whereas the meridional flow speed 
increases ~ i^q. These scalings change if i^o becomes so 
large that the viscous stress becomes a force comparable 
to the Coriolis force in the meridional flow equation, or 
if Kq becomes so small that advection of entropy by the 
meridional flow dominates over the diffusion. 



Solar differential rotation and meridional flow 



11 





0.0 V-: ........ 

0.0 0.2 0.4 0.6 0.8 1.0 
r/R 



c 
a. 



450 


(e) 




==^ 


— ; 




\ 


^~~_- 


400 


^ 


t^ 




- 




~ ~ — _ - 


^ — - ■ 






\^^^ 


— ^ ^_^^ 




350 






^^-^ 


^ 


300 











1 


— ^^i^^^i' . . 1 . . .^ 1 ■ 


, 1 , , , 




---^'>y^,^/' 


(c) 


0.8 
0.6 


^|i| 


}()<■ 


0.4 


v^ 


/' K- 


0.2 


- / „..■•'■■ Y 


// V 


0.0 


/,•-; ,,,,,,,, ,li 


/,,,,' 



0.0 0.2 0.4 0.6 0.8 1.0 
r/R 



c 

a 



450 


(0 ^^_— ^^=.== 




¥^JTrT^^j^^^^ : 


4UU 


V\ ' 




V""-C_:i^ — ' ^ ■ 




\ '"'' 


350 


^^^^^^^^^ 


300 





0.60 



1.00 



0.60 0.70 



0.80 
r/R 



0.90 1.00 



0.60 0.70 



0.80 
r/R 



0.90 1.00 



Fig. 6. — Solutions with a nonadiabatic convection zone, (a, d) Case 8 with r^^^, = 0.8, (b, e) case 9 with rgut, = 0.825, and (c, f) case 
10 with ^5^^, = 0.85. A superadiabatic convection zone can overcompensate for the effect of the subadiabatic tachochne, leading again to 
cylindrical differential rotation in the convection zone. However, if more than the lower third of the convection zone is weakly subadiabatic, 
this is sufficient to compensate for the upper superadiabatic part of the convection zone. 



4. SOLUTIONS WITH NONADIABATIC CONVECTION 
ZONES 

The models discussed so far assumed an adiabatic 
convection zone. This allowed us to separate the ef- 
fects of a subadiabatic tachocline from processes origi- 
nating within the convection zone such as the rotational 
anisotropy of convection, which are not considered in this 
model. 

Since the meridional flow generates additional entropy 
variations in a nonadiabatic convection zone, a consid- 
eration of a more realistic stratification within the con- 
vection zone is crucial, especially since a superadiabatic 
stratification generates entropy variations of opposite 
sign. Since we are using an axisymmetric mean field ap- 
proach, incorporating a superadiabatic stratification has 
the potential to introduce convective instability, which 
cannot be addressed appropriately in this model. We 
therefore have to make sure that with the turbulent vis- 
cosity and thermal diffusivity values we use, the Rayleigh 
number remains subcritical for convection. For the su- 
peradiabaticity profile defined in Eq. H25(l together with 
a value oi S = 3xl0~^at the top surface we have 
to increase the thermal diffusivity in the upper part of 
the convection zone to values of 5 x 10^ m^ s~^ to avoid 
convective instability. We suppress here the convective 



instability only by increasing Kt, since a change of vt 
would also alter the turbulent angular momentum trans- 
port and therefore change the differential rotation and 
meridional flow pattern in a way that makes a compari- 
son with the models discussed before difficult. Since we 
only change Kt above rsub the entropy perturbation that 
originates from the tachocline is only marginally affected. 

Nonlocal mixing-length models as discussed by 

Skale v fc Stixl l)1991il show typically below r = 0.75 to 
0.8 i?o a weakly subadiabatic convection zone with val- 
ues of J « —5 X 10^^. At r — 0.95 i?0 the superadia- 
baticity reaches values around 5 w 10~^ and increases 
then strongly toward the surface. We cannot represent 
the large surface values for the reasons described above; 
however, as we show later, they are also not that relevant 
for this problem. 

In Fig. 1^1 we shows models similar to case 1 but in- 
cluding a nonadiabatic convection zone as defined by Eq. 
(|25|) . We have varied here the parameter r^uh, which de- 
termines where the convection zone turns weakly suba- 
diabatic. In case 8 with rg^h = 0.8 Rq (panel a) and 
d)), the entropy perturbation arising from the superadi- 
abatic part of the convection zone is strong enough to 
overpower the effect of the subadiabatic tachocline to a 
large extent. Below 45° latitude the differential rota- 
tion is again close to the Taylor-Proudman state with 
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Fig. 7. — Model for values of r-gub between those in cases 9 and 10, with a nonvanishing inward angular momentum flux at the surface 
leading to a surface shear layer and an increased meridional return flow at the surface. In (a) we have overplotted dashed lines with a 25° 
inclination with respect to the rotation axis, which is the observed inclination of the Q contours between 15° and 55° latitude. 



cylindrical contour lines; at higher latitudes the disk like 
rotation contours are preserved, although the amplitude 
of the differential rotation is significantly reduced. Using 
a value of rgub — 0.85 i?© (case 10), the contribution of 
the subadiabatic part of the convection zone dominates 
over the superadiabatic top part, and a solution close to 
that of case 1 is retained (with slightly larger differential 
rotation). Comparing case 9 and case 10 shows that for 
the superadiabaticity profile we have chosen, the tran- 
sition takes place somewhere between j'sub — 0.825 and 
0.85 i?0. 

We found in our model that the role of a nonadiabatic 
convection zone depends on the extent of the weakly sub- 
adiabatic layer at the base. A transition in the solution 
takes place if roughly the lower 40% — 50% of the convec- 
tion zone is subadiabatic, even though the subadiabatic- 
ity in the lower part is very weak. Before we discuss 
whether such a solution is feasible at all (see next para- 
graph) , we investigate further how sensitively this result 
depends on the assumption of the superadiabaticity pro- 
file within the convection zone (e.g. we did not consider 
the strongly superadiabatic values close to the top of the 
convection zone). 

Eq. H20I) expresses clearly that the large values of the 
superadiabaticity near the surface do not matter that 
much, since the product g^S is relevant. The entropy 
perturbation at a given latitude in the convection zone 
results from an efficient diffusive exchange process of the 
quantity qoTq-si in radius, which means that the lower 
part of the convection zone contributes significantly, even 
though the value oi \S\ is very small there. Consequently, 
ii 6 ^ Q^ , which resembles the order of magnitude vari- 
ation in mixing- length models, then the contribution of 
the lower half of the convection zone is roughly of mag- 
nitude equal to the upper half. Therefore, if the lower 
half of the convection zone is weakly subadiabatic, ap- 
proximately — 10~^ to —10^^, this would dominate over 
the contribution of the superadiabatic upper part of the 
convection zone. 

Nonlocal mixing-length models do not provide such 
large s ubadia batic fractions of the convection zone. 
ISpruili ^^1997^ suggested a highly nonlocal convection 
model, in which the convection is driven in a narrow su- 



peradiabatic surface layer and downflows penetrate with 
a small amount of mixing all the way to the base of the 
convection zone. In this case the major fraction of the 
convection zone would be subadiabatic (because of radia- 
tive heating of the broad upffow in the lower convection 
zone) , but it is unclear to what extent this model applies 
to solar convection. 

In Fig. Ul we show a solution with a value of r^uh 
between the values used for case 9 and 10 (rgub — 
0.8375 Rq) and a different boundary condition for the 
angular velocity at the surface. Instead of setting the 
turbulent angular momentum transport to zero as ex- 
pressed in Eq. I|33() . we use a radially inward transport 
by setting A — —9 in a thin surface layer. Because of the 
stress-free boundary condition for Q this requires a neg- 
ative radial gradient of f2 near the surface. Th is param- 
eterization reflects the idea bv iGilman fc Fouk al (197!|) 
that the rotational influence on supergranulation leads 
to an outwardly decreasing Q. 

Fig. panel a) shows the differential rotation contours 
for this case. Over plotted are lines indicating a 25° in- 
clination angle with resp ect to the axis of rotat ion as 
found by hehoseismology i)Schou et al. 111998112001 . This 
feature is re produced very we l l from about 15° to about 
60° latitude. IGilman fc Howd l)2003|) presented an expla- 
nation for this phenomenon based on the influence of a 
one-cell meridional flow on the differential rotation, as- 
suming that the differential rotation would be constant 
in radius in the absence of the meridional flow. In our 
model it is difficult to identify the solution that would 
correspond to a solution in the absence of the meridional 
flow since the meridional flow is an integral part of the 
solution. Neglecting any effect of the meridional flow 
would provide a solution in which the Reynolds stress 
is balanced by diffusion, which in this particular case 
would yield contour lines with a 75° angle to the axis 
of rotation. In that sense the inclination is a result of 
the meridional flow but also including in a more compli- 
cated way the influence on the entropy profile within the 
convection zone. 

In panel c) we show the meridional flow speed at 45° 
latitude. A comparison with Fig. 0] shows an increase 
in the surface flow speed, which is a direct result of the 
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changed boundary condition for the angular velocity. 

We present this solution here in order to show that 
within the framework of this model it is possible to obtain 
solutions that are very close to the observed pattern by 
making assumptions about the superadiabaticity that go 
beyond the predictions of mixing-length theory (subadi- 
abatic part of convection zone extends up to 0.8375 Rq), 
but which are physically feasible if the degree of nonlo- 
cality of the convection is large enough. 

5. SUMMARY 
The main results of our model are as follows: 

1. The proflle of the differential rotation is deter- 
mined mainly by the profile of the entropy pertur- 
bation originating in the subadiabatic tachocline 
and spreading into the convection zone because of 
thermal conductivity. 

2. The profile of the differential rotation is rather in- 
sensitive to the parameterization of the Reynolds 
stress as long as there is a sufficiently large equator- 
ward angular momentum flux; however, the magni- 
tude of differential rotation changes with different 
assumptions. 

3. The parameterization of the Reynolds stress 
strongly inffuences the meridional ffow (compare 
cases 1, 2, 6, and 7). 

4. If the lower half of the convection zone is weakly 
subadiabatic, the solar differential rotation can be 
explained through entropy perturbations arising 
from the nonadiabatic stratification. If the sub- 
adiabatic region has a smaller extent, additional 
effects such as anisotropic heat transport are re- 
quired. 

5. For angular momentum transport almost aligned 
with the axis of rotation and angular momentum 
transport in latitude only, we find a dominant coun- 
terclockwise meridional ffow cell (equatorward at 
the base of and poleward at the upper layers of 
of convection zone) as a robust result. Our model 
shows also the tendency of a weaker reverse cell 
at higher latitudes, which has been observed by 
Ulaber et al. (2002). A dominant counterclockwise 
meridional flow cell is favorable for flux-transport 
dynamo models, in which the equatorward merid- 
ional flow at the base of the convection zone ensures 
the equatorward propagation of magnetic activity 
through the solar cycle. 

6. The meridional flow shows very little penetration 
beneath the base of the convection zone. In ad- 
dition to the constraints imposed by the subadia- 
batic stratification in the radiative interior, angular 
momentum conservation does not allow for a sig- 
nificant meridional flow there, since the associated 
angular momentum transport would lead to signifl- 
cant changes of differential rotation unless opposed 
by a strong Reynolds stress (which is highly un- 
likely in the radiative interior). 



6. IMPLICATION FOR SOLAR DIFFERENTIAL ROTATION 

We presented in this paper a simplified model for the 
solar differential rotation and meridional flow. This 
model parametrizes important convective scale processes 
such as the turbulent angular momentum transport and 
turbulent diffusivities and does exclude processes such 
as rotational anisotr opy of the convective energy flu x as 
discussed in d e tail bv lKitchatinov fc Riidigcr (1995) and 
iKiiker fc Stixl l)200ll) . This model is therefore not in- 
tended to be a complete differential rotation model for 
the solar convection zone, but rather a model to evaluate 
the importance of effects resulting from a nonadiabatic 
stratification in tachocline and convection zone. 

Even though we did not use here exact ly the formula- 
tion of the A-effect as a dopted bvlKitchati nov fc R.udigeil 
iITqqI . IRiidiger et aP iITqqI . and .Kiiker fc Stixl (120011^ 
for the solar case, our results are in general agreement 
with their earlier work. For a magnitude of the A-effect 
(the parameter Ag in our model) of order unity, we get an 
amplitude of the differential rotation comparable to solar 
values. Differences occur in the profile of the differential 
rotation (inclination of isolines with respect to the axis 
of rotation) because of the different physics considered 
in the entropy equation of our model (which is the main 
focus of this work). Larger, but explainable differences, 
exist in the meridio nal flow patterns obta ined. Although 
the investigation of iRiidieer et alJ l)1998j) shows a coun- 
terclockwise flow similar to our results in a model includ- 
ing a significant amount of viscosity and no latitudinal 
entropy variation, models including a latitudinal entropy 
gradient and a A-effect formulation taking into account 
the variation of the Coriolis number within the convec- 
tion zone typically yield a clockwise flow cell close to 
the s urface (Kitchatinov fc Riidiger 1995; Kiiker fc Sti3 
1200 1|) . As explained by the authors this is because of a 
radially outward transport of angular momentum in the 
upper layers, which is not parameterized in our model. 
For most of our cases we assume a vanishing angular mo- 
mentum flux at the surface; in the case shown in Fig. d 
we use a radially inward transp ort of angular momentum 
(as suggested by the work of iGilman fc FoukaJI 1)19791) 
and also by observations, which clearly show a decrease 
of rotation rate with radius in the outer layers), which has 
the opposite effect of enhancing the poleward meri diona l 
surface flow. The investigation of iKiiker fc Stixl l|200J) 
also indicates that the meridional flow is much more sen- 
sitive to details of the model used than the differential 
rotation (see Fig. 2 therein, which addresses the influ- 
ence of the mixing- length parameter), in agreement with 
the flndings in this paper. 

The models with adiabatic convection zones, which are 
shown in Fig. El to allow evaluating of the impor- 
tance of the tachocline effects for models such as that 
of Kiiker fc Stix ( 2001) that only consider the convec- 
tion zone. To summarize our results, the entropy per- 
turbation generated within a subadiabatic tachocline is 
strong enough to avoid the Taylor-Proudman state above 
30° latitude and is therefore as important as other ef- 
fects such as rotational a nisotropy of the convect i ve en- 
ergy transport. However. Kitch atinov fc R.iidigeil J1995) 
found in their investigation that the effect of a subadi- 
abatic tachocline is rather small compared with the ef- 
fect of anisotropic energy transport. This apparent con- 
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tradiction could be caused either by the diflFerent steep- 
ness of the shear profile o f the differential rotation in the 
tachocline in the model of lKitchatinov fc Riidieeil l)1995D 
(see, e.g., Fig. 1 and 2 therein), since the magnitude of 
the entropy perturbation depends strongly on the value 
of dfl/dz, or by an insufhcient overlap between thermal 
conductivity profile and subad iabatic overshoot region 
IjKitchatinov &: Riidieeil l)1995 f) computed the turbulent 
diffusivity through a local mixing- length relation). 

Since the value of di^/dz within the subadiabatic 
tachocline is significantly larger than the value of dfl/dz 
in the convection zone, the magnitude of the expected 
entropy perturbation in the tachocline also exceeds the 
magnitude required to balance a solar like differential 
rotation within the convection zone. If therefore only a 
fraction of the entropy perturbation generated within the 
subadiabatic tachocline spreads into the convection zone, 
this provides a significant contribution. This spread 
depends mainly on the overlap between the subadia- 
batic tachocline and the region that is mixed by con- 
vection. In our model we typically use a convective dif- 
fusivity of 5 X lO^m^s"^ at r — rtran = 0.725 i?© and 
lower values below, which introduced enough coupling 
between the subadiabatic region and the convection zone. 
Since the observed tachocline spreads about one-third 
of the way into the region above r = 0.713 i?©, which 
is very close to adiabatic according to helioseismology 
IjCharbonneau et al]ll999j) . it can be expected that there 
is a sufficient coupling between these two regions in the 
case of the Sun. 

The models with a nonadiabatic convection zone (Figs. 
El ITJ allow us to estimate under which conditions the en- 
tropy perturbation resulting from the nonadiabatic strat- 
ification is sufficient to explain the observed solar dif- 
ferential rotation without any additional effects such as 
anisotropic heat transport. The main problem is that a 
superadiabatic convection zone overcompensates the ef- 
fect of the subadiabatic tachocline to a large degree un- 
less the lower half of the convection zone is weakly suba- 
diabatic. Even though the absolute values of the supera- 
diabaticity in the convection zone are much lower than 
those in the tachocline, their contribution can be very 
large because of the much larger radial meridional flow 
velocity within the convection zone. We found that the 
effect of the superadiabatic part of the convection zone 
can be compensated by a weakly subadiabatic lower half 
of the convection zone, since the thermal inertia of the 
entropy perturbation is r^ Po^oSi and therefore even a 
small entropy perturbation at the base of the convection 
zone can contribute more than a large entropy perturba- 
tion in the upper part of the convection zone. Within the 
frame work of this model we cannot address the question 
of whether such a solution is feasible, since this would 
require more sophisticated nonlocal convection theory. 
Nonlocal mixing-length models typically predict a suba- 
diabatic stratification below r = 0.8 i?©; however, for a 
larger degree of nonlocality (stable downflows with only 
little entrainment and detrainment) a larger extent of 
this region would be possible. If additional effects such 



as anisotropic heat transport and nonlinear feedback of 
the stratification on the heat conductivity were included, 
it could be possible that a solar-like solution as shown in 
Fig. [3 is possible with a smaller extent of the subadia- 
batic part of the convection zone. 

A further test of the process proposed in this paper 
would be the inclusion of a subadiabatic tachocline in 
the full spherical shell convection simulations. Besides 
adding a subadiabatic overshoot region, this requires the 
inclusion of the tachocline shear layer either by resolving 
the relevant physical processes or by forcing the shear (as 
done in this model) through a lower boundary condition, 
since the entropy perturbation arises as a consequence of 
a subadiabatic shear layer. Since the entropy perturba- 
tion is also affected significantly by the meridional flow 
within a nonadiabatic convection zone, it is also crucial 
to get the meridional flow pattern right. Work in this 
direction is currently in progress. 

Since the meridional flow is of high interest for flux 
transport dynamos, here we discuss in more detail the 
findings of our model concerning the penetration of the 
return flow below the base of the convection zone. As 
mentioned above, our model provides a counterclockwise 
flow as a robust result, provided that there is an inwardly 
directed turbulent angular momentum flux. The return 
flow at the base of the convection zone is located in the 
region with the largest turbulent viscosity gradient, since 
there the divergence of the turbulent angular momentum 
flux (A-effect) would lead to an increase of rotation rate 
unless opposed by an equatorward transport of angu- 
lar momentum by the meridional flow (which tends to 
slow down the rotation rate). In that way the merid- 
ional flow is always tied to the presence of turbulent 
Reynolds stress. For the parameterization of the tur- 
bulent diffusivities we use in most of our models, the 
meridional flow speed at r = 0.71 Rq typically drops to 
less than 10% of the maximum return flow speed, which 
is reached at around r = 0.745 i?©. Below r = 0.71 i?© 
velocities on the order of a few cms^^ persist, how- 
ever our model tends to overestimate their amplitude 
since we use overshoot values for 6. Using radiative core 
values would decrease this amplitude further by several 
orders of magnitude. Therefore, a penetration of the 
meridional flow below the base of the convection zone 
as used in a few flux transport dynamo models seems 
very unreasonable, since it implies that there are very 
strong Reynolds stresses in the strongly subadiabatic ra- 
diative core. Since the constraint set by the angular mo- 
mentum conservation is additional to the constraint set 
by the subadiabaticity of the stratification, which was 
pointed out by Gilman & Micsch (2004), the limitation 
on penetration of meridional flow beneath the base of 
th e convection zone i s even more stringent than found 
bv lGilman fc MieschI l|200|) . 



The author wants to thank P. A. Gilman and M. S. 
Miesch for stimulating discussions about the solar differ- 
ential rotation problem. 
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